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(57) Abstract: The invention features a cali- 
bration method for diffuse optical measurements 
that corrects trans mittance measurements 
between a source and a detector for factors 
unrelated to sample properties. The calibration 
method is based on the same set of transmittance 
measurements that are subsequently corrected 
by the calibration and used in imaging and/or 
spectroscopy applications. The calibration 
method involves a forward calculation for each 
of multiple source-detector pairs based on an 
approximate model of die sample (310), and a 
minimization of an expression (330) that depends 
on the forward calculations and the transmittance 
measurements to determine self-consistent 
coupling coefficients for every source-detector 
pair. Once the coupling coefficients have been 
detennined (360), they can be used to correct 
the transmittance measurements. If desired, an 
inverse calculation (370) can be performed on 
the corrected sample measurements to determine 
spatial variations in the optical properties of 
the sample. If necessary, the calibrat tion can 
be repeated (380) and iteratively improved, 
whereby the optical properties determined by the 
inverse calculation (370) in an earlier iteration 
are used to improve the sample model for the 
forward calculation in a subsequent iteration. 
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CALIBRATION METHODS AND SYSTEMS FOR 
DIFFUSE OPTICAL TOMOGRAPHY AND SPECTROSCOPY 
Cross-Reference to Related Application 
This application claims the benefit of U.S. Provisional Application No. 
5 60/1 54,423, Tiled on September 17, 2000, which is incorporated herein by reference in its 
entirety. 

Field of the Invention 
The invention relates to the calibration of optical techniques for imaging and 
spectroscopy of, e.g., biological tissue. 

10 

Background of the Invention 
Recently there has been significant interest in using optical radiation for 
imaging within highly scattering media, such as biological tissue. Photons travel within 
the highly scattering media along a distribution of paths, of which very few are straight. 

1 5 Thus, directing light into the highly scattering media and subsequently detecting the 
diffuse light emitted from the media provides information about local variations in the 
scattering and absorption coefficients. Such information can identify, for example, an 
early breast or brain tumor, a small amount of bleeding in the brain, or an early aneurysm. 
Furthermore, for example, multiple wavelengths can be used to spectroscopically 

20 determine local tissue concentrations of oxy-hemoglobin (HbO) and deoxy-hemoglobin 
(Hb) in tissue, which may be in response to some stimulus, e.g., a drug. For a general 
description of such applications, see, e.g., A. Yodh and B. Chance in "Spectroscopy and 
Imaging with Diffusing Light," Physics Today, pp. 34-40 (March 1995). 

If the spatially varying optical properties of the highly scattering media are 

25 known, photon propagation within the media can be calculated numerically. The 

numerical calculation is simplified when scattering is much larger than absorption, in 
which case the photon propagation can be approximated as a diffusion process. This 
condition is typically satisfied in biological tissue in the spectral range of about 700 nm to 
900 nm. The numerical calculation gives the distribution of light inside the tissue, and is 

30 usually referred to as the "forward calculation." For a sample being imaged, however, the 
"inverse calculation" needs to be solved, i.e., deducing the sample's optical properties 
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from the diffuse light measurements. Numerical techniques for performing the inversion 
include singular value decomposition (SVD), simultaneous iterative reconstruction 
technique (SIRT), K-space diffraction tomography, and using an extended Kaman filter. 
For a general review of techniques for the forward and inverse calculations, see, e.g., S.R. 
5 Arridge in "Optical tomography in medical imaging," Inverse Problems, 15:R41-R93 
(1999). 

In diffuse optical tomography (DOT), multiple sources sequentially direct 
light into tissue at spatially separated locations, and for each such source, multiple 
detectors on the tissue measure the diffuse light emitted from the sample at spatially 

10 separated locations. For every source-detector pair, one measures the local transmittance 
of the diffuse light, i.e., the ratio of the diffuse radiance measured by the detector and the 
incident radiance from the source. The measurements provide the input information for 
the inverse calculation. However, the measurements can include various uncertainties 
caused by, for example, fluctuations in the source power, variations in the detector gain, 

1 5 and coupling variations at the source-tissue interface as well as the tissue-detector 
interface. 

To minimize the uncertainties, DOT systems are typically calibrated with 
initial measurements for a known sample, and the calibration is used to correct 
subsequent measurements for imaging an unknown sample. Unfortunately, coupling at 

20 the source-tissue interface and the tissue-detector interface can vary from measurement to 
measurement because of, for example, movement or perspiration of the patient, or 
movement of an optical fiber that forms part of a source or detector. Thus, the results of 
the inverse calculation can include systematic errors caused by measurement variations 
that are independent of the tissue optical properties of interest. 

25 The systematic errors can also limit absolute spectroscopic measurements of 

optical properties at a particular spatial location, e.g., the absolute, rather than relative, 
values of absorption and scattering. 

Summary of the Invention 
30 The invention features a calibration method for diffuse optical measurements 

that corrects transmittance measurements between a source and a detector for factors 
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unrelated to sample properties. For imaging applications, the corrected transmittance 
measurements can be subject to an inverse calculation to determine spatial variations in 
the optical properties of the sample, i.e., to "image" the sample. For spectroscopic 
applications, the corrected transmittance measurements can be used to determine absolute 

5 values for the optical properties of the sample in a particular spatial region at multiple 

wavelengths, e.g., to determine the absolute concentrations of oxy-hemoglobin (HbO) and 
deoxy-hemoglobin (Hb). The calibration method is based on the same set of 
transmittance measurements that are subsequently corrected by the calibration and used in 
imaging and/or spectroscopy applications. The accuracy of the subsequent results is thus 

10 not subject to uncertainties caused by a delay between calibration and sample 
measurements. 

The calibration method involves a forward calculation for each of multiple 
source-detector pairs based on an approximate model of the sample, and a minimization 
of an expression that depends on the forward calculations and the transmittance 

1 5 measurements to determine self-consistent coupling coefficients for every source-detector 
pair. Once the coupling coefficients have been determined, they can be used to correct 
the transmittance measurements. If desired, an inverse calculation can be performed on 
the corrected sample measurements to determine spatial variations in the optical 
properties of the sample. If necessary, the calibration can be repeated and iteratively 

20 improved, whereby the optical properties determined by the inverse calculation in an 
earlier iteration are used to improve the sample model for the forward calculation in a 
subsequent iteration. 

In general, in one aspect, the invention features a system for making optical 
measurements. The system includes at least two optical sources which during operation 

25 couple optical radiation into a sample at spatially separated locations and at least two 
optical detectors positioned to receive optical radiation emitted from the sample at 
spatially separated locations in response to the optical radiation from the sources, and an 
analyzer. 

The signal g(ij) produced by the fit detector in response to the optical 
30 radiation from the / th source can be expressed as g(ij) = S 1 D- '/(/,./), where f(ij) 
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depends only on the properties of the sample, 5 s is a coupling coefficient for the i* 
source, and Z)j is a coupling coefficient the detector. During operation, the analyzer 



calculates the value of the product SlD k for at least one of the source-detector pairs based 
on the signals produced by the detectors and simulated values oT/OJ) corresponding to a 
5 model of the optical properties of the sample. The sources, for example, can provide 
continuous-wave radiation, in which case g(ij),f(ij), and D) are all real-valued. 
Alternatively, the sources can provide modulated CW radiation, or short temporal pulses 
of optical radiation (e.g., less than about 1 to 100 ps). In these latter cases, g(ij),f(ij), S 1 , 
and Di can be complex, or more generally they can be vectors representing multiple 
1 0 values (e.g., transmittance and temporal delay). 



detector pair based on the detector signals and the simulated values of f(ij). The analyzer 
can further calculate experimental values of f(ij) based on the calculated values of S*D1 

1 5 and the signals g(ij) using the expression g(ij) = 5'Z) y /(/,y), and then perform an 
inverse calculation on the experimental values for f(ij) to determine spatial variations in 
at least one optical property of the sample. The optical property or properties can be the 
absorption coefficient, the reduced scattering coefficient, or both. The analyzer can 
modify the model of the sample based on the determined spatial variations, and then 

20 repeat the calculation of the values of the product for every source-detector pair 
using the modified model. The initial model can correspond to the sample being 
homogeneous. 

The analyzer can simulate values of f(ij) according to the expression: 



Embodiments of the system can also include any of the following features. 
The analyzer can calculate the value of the product &DJ for every source- 



25 




where jr tf | is the distance between the I th source and the j* detector. This expression 



corresponds to an infinite, homogenous medium. 
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The analyzer can calculate the value of the product SlD k by minimizing the 

expression: 



F(S'D> ) = ±t{ L ° s i D k : 1 } fOJ) - g(U)j . 



where 



L(i,j,k,l) = A?Ai' 



,k A'.v ^ g{ij) 



A J> = 



and Ns is the number of sources and No is the number of detectors. 
10 Furthermore, when the model corresponds to the sample being homogeneous, 

the analyzer can calculate at least one of the absorption coefficient \± a and the reduced 
scattering coefficient \x.' s by minimizing F{s'D k ) with respect to the product S'Z)k and 
the at least one of the absorption and scattering coefficients. F(s'D k ) implicitly depends 
on \i a and \x! s through f(ij). Similarly, the analyzer can calculate both of the absorption 

1 5 and scattering coefficients by minimizing f(s'd") with respect to the product S^D^ and 
the absorption and scattering coefficients. 

The analyzer can also calculate the product S™D n for every source-detector 
pair by minimizing F(s m D n ). Alternatively, the analyzer can calculates the product SPDi 

for every remaining source-detector pair according to S'D J = ^^/'*'^ . 

20 In another aspect, the invention features a method for calibrating an optical 

measurement system. The optical measurement system includes at least two optical 
sources and at least two optical detectors. The sources couple optical radiation into a 
sample at spatially separated locations and the detectors are positioned to receive optical 

-5- 



WO 01/19241 



PCT/US00/25467 



radiation emitted from the sample at spatially separated locations and generate signals in 
response to the optical radiation from the sources. 

The methods includes the steps of: providing the signals generated by the 
detectors, wherein the signal g(ij) generated by the 7 th detector in response to the optical 
5 radiation from the f th source can be expressed as g(i 9 j) = S'D J f(i,j), where f(ij) 
depends only on the properties of the sample, S* is a coupling coefficient for the i* 
source, and Z>J is a coupling coefficient for the detector; and calculating the value of 
the product S^D^ for at least one of the source-detector pairs based on the signals 
generated by the detectors and simulated values of f(ij) corresponding to a model of the 
1 0 optical properties of the sample. 

The sources in the optical measurement system, for example, can provide 
continuous-wave radiation, in which case g(ij),f(ij), and DJ are all real-valued. 
Alternatively, the sources can provide modulated C W radiation, or short temporal pulses 
of optical radiation (e.g., less than about 1 to 100 ps). In these latter cases, g(ij),f(ij) y 
15 and BS can be complex, or more generally they can be vectors representing multiple 
values (e.g., transmittance and temporal delay). 

Embodiments of the calibration method can also include any of the following 

features. 

The calibration method can calculate the value of the product S*D3 for every 
20 source-detector pair based on the detector signals and the simulated values of f(ij). The 
method can further include calculating experimental values of f(ij) based on the 

calculated values of SiCU and the signals g(ij) using the expression g(ij) = S'D J f(i,j), 
and then performing an inverse calculation on the experimental values for f(ij) to 
determine spatial variations in at least one optical property of the sample. The optical 
25 property or properties can be the absorption coefficient, the reduced scattering coefficient, 
or both. Furthermore, the method can include modifying the model of the sample based 
on the determined spatial variations, and then repeating the calculation of the values of 
the product SlDJ for every source-detector pair using the modified model. The initial 
model can correspond to the sample being homogeneous. 
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The simulated values of ffiJJ can be calculated according to the expression: 



A**j)<* — n ex P 
4 T</I 



-(3n> fl )i|ry| 



where |r^ | is the distance between the i* source and the j th detector. This expression 
corresponds to an infinite, homogenous medium. 

The method can calculate the value of the product S^D^ by minimizing the 

expression: 



where 

10 L(i,j,kJ) = A?Aj 





gO,j) 












gOJ) 







and Ns is the number of sources and N D is the number of detectors. 
1 5 Furthermore, when the model corresponds to the sample being homogeneous, 

the method can include calculating at least one of the absorption coefficient ]i a and the 

reduced scattering coefficient ii' s by minimizing F{s'D k ) with respect to the product 

£'Z) k and the at least one of the absorption and scattering coefficients. f(s'd") 

implicitly depends on \x a and \l s through f(ij). Similarly, the method can include 

20 calculating both of the absorption and scattering coefficients by minimizing F(s'D k ) 

with respect to the product S'D^ and the absorption and scattering coefficients. 



WO 01/19241 



PCT/US00/25467 



The method can also include calculating the product 5 m Z) n for every source- 
detector pair by minimizing F(s m D n ). Alternatively, the method can calculates the 

product SiDJ for every remaining source-detector pair according to S i D j = MLl^D. 

S'D k ' 

In another aspect, the invention features a computer-readable medium which 
5 causes a processor to perform the steps of any of the embodiments of the calibration 
method described above. 

Unless otherwise defined, all technical and scientific terms used herein have 
the same meaning as commonly understood by one of ordinary skill in the art to which 
this invention belongs. Although methods and materials similar or equivalent to those 
10 described herein can be used in the practice or testing of the present invention, the 

suitable methods and materials are described below. All publications, patent applications, 
patents, and other references mentioned herein are incorporated by reference in their 
entirety. In case of conflict, the present specification, including definitions, will control. 
In addition, the materials, methods, and examples are illustrative only and not intended to 
15 be limiting. 

Embodiments of the invention include many advantages. For example, the 
new calibration methods eliminate errors caused by using separate data sets to calibrate 
and image samples. Where the sample is human tissue, such errors can result from 
patient or fiber movement, or patient perspiration. Furthermore, the calibration method 
20 obviates the need to minimize fluctuations in the source outputs and detector gains. 
Diffuse optical measurement systems employing the new calibration methods can 
therefore provide more accurate and robust results. 

Other features and advantages of the invention will be apparent from the 
following detailed description, and from the claims. 

25 

Brief Description of the Drawings 
FIG. 1 is a schematic diagram of a diffuse optical tomography (DOT) system. 
FIG. 2 is a schematic diagram of a source-detector pair for the system of FIG. 

1. 
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FIG. 3 is a flow chart summarizing the steps of one embodiment of the 
calibration method described herein. 

FIG. 4 is a schematic diagram of electronic hardware for the DOT system of 

Example 1. 

5 FIG. 5 is a schematic diagram of the optical arrangement of the sources and 

detectors of Example 2. 

FIGS. 6a-6b are graphs of optical coefficients determined in Example 2 based 
on experimental data. 

FIGS. 7a-7d are graphs of optical coefficients determined in Example 3 based 
1 0 on simulated data. 

FIGS. 8a-8b are reconstructed DOT images determine with (FIG. 8b) and 
without (FIG. 8a) the calibration method described herein for the simulation of Example 
4. 

FIG. 9 is a schematic diagram of the optical arrangement of the sources and 
1 5 detectors of Example 5 . 

FIGS. 10a- lOd are reconstructed DOT images of experimental data from 

Example 5. 

FIGS. 1 la-b are graphs of peak image values derived from the images of 
FIGS. lOa-lOd. 

20 FIG. 12 is an additional reconstructed image of the experimental data from 

Example 5. 

Detailed Description 

General System 

The invention features a calibration method for a diffuse optical measurement 
25 system, e.g., a diffuse optical tomography (DOT) system. A schematic diagram of a DOT 
system is shown in FIG. 1. The system 100 includes an array of spatially separated light 
sources 1 10 and spatially separated detectors 120. During use, the array of sources and 
detectors is positioned over a sample 1 50 to be imaged, e.g., a patient's head or breast. A 
controller 130 connected to light sources 110 sequentially triggers them to couple light 
30 into sample 150, which is a highly scattering media that causes the light to become 
diffuse within the sample. For each sequentially triggered source, each detector 120 

-9- 
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measures the light that reaches it through sample 150. Controller 130 is also connected to 
detectors 120 and selectively channels the signals from the detectors. An analyzer 140 is 
connected to controller 130 and analyzes the signals measured by detectors 120. 

The signal g(ij) measured by the / h detector in response to light coupled into 
5 the sample by the fa source can be expressed as: 

g{i,j) = S'Djf(i,j) (1) 

where f(ij) is the transmittance of the sample from the fa source to the ft* detector, and 

1 0 S* and DJ are coupling coefficients for the fa source to the /h detector, respectively. The 
transmittance f(ij) depends only on the optical properties, e.g., the spatially varying 
absorption and scattering coefficients, and can be numerically calculated by a forward 
calculation if the optical properties are known. Conversely, if the transmittance^^ can 
be calculated from the measured signals g(ij), an inverse calculation can be performed on 

1 5 f(UJ to yield the optical properties of the sample and reveal, e.g., the presence of an 

object hidden in the sample. The source coupling coefficient S* includes all of the factors 
associated with coupling light generated from the fa source into the sample. The detector 
coupling coefficient DJ includes all of the factors associated with coupling light out of the 
sample to generate an output signal at the fa detector. 

20 For example, in one embodiment shown in FIG. 2, each source 120, e.g., the 

fa source, includes a diode laser 200 for producing the optical radiation, an optical fiber 
210, and a lens 220 for coupling the optical radiation into fiber 210, which includes an 
end 215 adjacent sample 150 for directing the optical radiation into the sample. Each 
detector 130, e.g., the / h detector, includes an optical fiber 250 having an end 255 

25 adjacent sample 150 for receiving the optical radiation emitted from sample 150, a 

photodetector 260 for measuring the intensity of the optical radiation received by fiber 
250, and an amplifier 270 for amplifying the output of photodetector 260 to give the 

measured signal g(ij). In this embodiment, the source coupling coefficient & is the 
product of the fluence P L produced by diode laser 200, the coupling coefficient C LF of the 
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lens 220 to the fiber 210, the transmission coefficient TLf of fiber 210, and the coupling 
coefficient Cpj at the interface of fiber 2 1 0 and sample 1 50. Similarly, the detector 
coupling coefficient Di is the product of the coupling coefficient C T f at the interface of 
fiber 250 and sample 150, the transmission coefficient T LF of fiber 260, and the coupling 
5 coefficient of the sample fluence produced by the diode laser P L , the coupling coefficient 
Cfd between fiber 260 and photodetector 260, the efficiency y PE of photodetector 260, 
and the gain A of amplifier 270. 

In other embodiments, the light source can include a laser other than a diode 
laser, e.g., an ultrafast laser, or instead it can include an incoherent source. Also, the 

10 sources can include a common light source that selectively couples light into one of 
multiple fibers that deliver the Jight to spatially separated locations on the sample. 
Alternatively, the sources need not include optical fibers at all, for example, the lasers 
themselves can be positioned adjacent the sample or can include beam delivery optics to 
direct the light to the sample through free space. Furthermore, the light sources can 

1 5 provide light at multiple wavelengths by including, e.g., multiple diode lasers. 

To analyze the measured values gfij), a calibration is required to convert the 
measured values for g(i J) into the transmittance values f(i J). According to Equation (1), 
this requires a calibration for the value of SPDS for every source-detector pair. As will be 
described below, analyzer 140 determines self-consistent values for S*D1 based on the set 

20 of measured values g(ij) and the results of a numerical calculation for f(ij) corresponding 
to an approximate model of the sample: Once the calibration is determined, the same set 
of measured values g(ij) can be used to calculate f(ij) according to Equation (1), from 
which the analyzer can perform the inverse calculation. Because the same set of 
measured values are used for the calibration and the inverse calculation, the optical 

25 properties determined by the analyzer do not include systematic errors caused by 
fluctuations in the source and detector coupling coefficients. 

In the description that follows, it is assumed that the sources provide CW 
optical radiation and the detectors measure the intensity of the optical radiation, in which 
case g(ij),f(ij) y tf, and £U are all real-valued. However, the calibration techniques 

30 described herein can also be applied to other diffuse optical measurement techniques in 

- 11 - 
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10 



20 



which the sources do not provide CW radiation. For example, in some techniques, the 
amplitude of the optical radiation provided by the source is modulated to create photon 
density waves in the sample, and the detectors are configured to measure the amplitude 
and phase of the photon density waves after propagation through the sample. In this case 
the values for g(ij),f(ij), S\ and DJ can be complex. For a general reference on DOT 
with modulated optical radiation see, e.g., M.A. O'Leary et ah, Phys. Rev. Lett. 69. ? 
(1992). Furthermore, in other techniques, each source provides a temporally cohc : 
light pulse, e.g., a picosecond pulse, and the detectors are time-gated to measure the 
temporal delay of the diffuse light pulse in addition to its intensity. For a general 
reference on such time-domain DOT techniques see, e.g., M.S. Patterson et ah, Appl Opt. 
28, 2331, (1989), and S.R. Arridge in "Optical tomography in medical imaging" ibid.. 



Calibration Method 

Assuming N s sources and detectors and following Equation (1 ), the 
1 5 coupling coefficient of the 7 th source and 7 th detector can be expressed as: 



s' = J-y gtti> (2 ), 



D > = J_y 80, J) (3) 



s . D j = giU) (4). 



Equation (2) is the average of all measurements made with source /' and Equation (3) is 
the average of all measurements made with detector j. Substituting Equations (3) and (4) 
25 into Equation (2) we obtain a relationship between the source coupling coefficient 5* and 

the detector coupling coefficient D\ corresponding to the /* source and the k th detector, 
respectively: 
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and likewise between 23) and 5', corresponding to the source and detector coupling 
5 coefficients for the ft 1 detector and the l tn source, respectively: 



D j = *LlY gOJ) = AS_ (6). 



10 



Thus, 



, _ L(i,j,k,l) 
S'D" 



(7), 



where 



15 L(iJ 9 k 9 l) = A' s k Af (8). 

Note that L(iJ,k t l) is a function only of the number of sources and detectors, N s and Nj, 
the measurements, g(ij), and the optical properties of the medium as reflected through 
fjj). In other words, L(iJ,k,l) is not dependent on the coupling coefficients. 
20 Combining Equations (2), (3), and (7), the coupling coefficient product SlD^ 

must minimize the following expression to be consistent with each other and with the 
measurements g(ij): 

F(S'D> ) = ££(MUgD. fiiJ) _ g(/ , y) Y ( 9 ). 

25 
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Assuming that^v) can be calculated from an approximate model, the value for S*D k can 
be calculated by minimizing Equation (9) with respect to S^D^. The minimization can be 
repeated for every source-detector pair, or alternatively, once the value for S^D^ is 
determined for the l-k source-detector pair, the other coupling coefficient products can be 
5 determined from Equation (7). Although the calibration is non-linear, with no unique 
solution for the individual coupling coefficients 5* and £U, an exact solution exists for the 
coupling coefficient product of every source-detector pair, which is all that is needed for 
the calibration. The minimization of Equation (9) can be performed using standard 
techniques known in the art, see, e.g., W.H. Press et al. in Numerical Recipes in C: The 
10 Art of Scientific Computing (Cambridge U. Press, New York, 1988). 

To perform the minimization in Equation (9), analyzer 140 makes a numerical 
forward calculation for based on an approximate model of the sample. An initial 
model can be that the sample is homogenous with a constant absorption coefficient \x a 

and a constant reduced scattering coefficient \x! s . Once the calibration is performed, an 
1 5 inverse calculation can provide perturbative corrections to the model to show spatial 
variations in the optical properties of the sample. Numerical techniques for the forward 
and inverse calculations are known in the art and will be briefly described in the next 
section. However, for the simple case of an infinite homogeneous sample, the model 
forward calculation simplifies to the following expression: 

20 

/(^^^^^exp^S^^^^I^I (10), 

where Jr,yj is the distance between the I th source and the detector on the sample surface. 

If the background optical properties are not known, then an iterative procedure 
25 can be used to find \i a , \x! s , and the coupling coefficients. The procedure involves 

estimating \i a and \i' s to calculate^), and then minimizing F^D*) in Equation (9) to 
find the coupling coefficient SfD^. The procedure is repeated with different values of \x a 
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and yi' s in order to minimize F(ja fl ,jj^ ySfD^) with respect to all three parameters. For 

spectroscopic applications where constant values of \x a and are expected over the 

spatial extend of the measurement, the technique provides an accurate determination of 

the absolute values of \x a and \x s . 

5 Once all of the coupling coefficients are determined based on the initial 

model calculation, experimental values for J[i J) are calculated from Equation (1) based on 
the determined coupling coefficients and the measurements g(ij). The analyzer then 
performs an inverse calculation on the experimental values for J(iJ) to determine 
perturbations to the homogeneous model for the sample. If necessary, the calibration can 
10 be repeated for an improved model of the sample based on the results of the inverse 
calculation. In turn, the inverse calculation can be repeated for experimental values 
calculated from Equation (1) using the revised calibration. This iterative process can be 
repeated until the results for the spatially varying optical properties of the sample begin to 
converge. 

1 5 The steps performed by the analyzer to carry out the calibration and analysis 

of measurements g{ij) for imaging applications are summarized by the flow chart in FIG. 
3. 

In step 300, the analyzer receives measured signals g(ij) from the detectors. 
In step 3 10, an initial model for the sample is input into the analyzer, for 
20 example, the initial model may treat the sample as being homogeneous. In this case, 

values for \x a and \x s are estimated if they are otherwise unknown. 

In step 320, the analyzer calculates values iorfJJ) based on the sample model 
and a forward calculation. For example, if the sample is modeled to be homogeneous and 
infinite, Equation (10) can be used. 
25 In step 330, the coupling coefficient SfD^ is determined by minimizing F(5> l/) 

with respect to SfD* in Equation (9), the other parameters in Equation (9) being specified 
by the measurements for g(ij) and the values fox fljj) from the model forward calculation 
in step 320. 
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In step 340, steps 3 10-330 are repeated as necessary with additional estimates 
for ii a and n' 5 > until values are found for SfD*, \x a , and \x! s that optimally minimize 

In step 350, the values for all of the remaining coupling coefficients are 
5 determined by either: 1) calculating using Equation (7) and the value for sfrf 

determined in step 340; or 2) replacing the argument SfD* for F in Equation (9) with the 
coupling coefficients SD? corresponding to each of the remaining source-detector and 
determining the value SD? that minimizes F. 

In step 360, the calibration coefficients Si? determined in steps 340 and 350, 
1 0 and the measurements gfJJ) from step 300 are used to calculate experimental values for 
fijj) based on Equation (1). 

In step 370, an inverse calculation is performed on the experimental values for 
J{iJ) calculated in step 360 to determine spatial variations in the optical properties of the 
sample, e.g., an object hidden with a highly scattering sample. 
15 1° ste P 380, if necessary, the model for the sample estimated in step 3 10 is 

revised based on the results from step 370, then steps 320, 330, and 350 are repeated, one 
time, to recalculate the calibration coefficients Stf based on the revised sample model. 

In step 390, if necessary, steps 360-370 are repeated, one time, using the 
recalculated calibration coefficients from step 380 to improve the determination of the 
20 spatial variations in the optical properties of the sample. Steps 380-390 can be iteratively 
repeated as necessary until the spatially varying optical properties determined in step 390 
converge to within a desired accuracy. 

For spectroscopic applications in which the optical properties of the sample 
are expected to be homogeneous, the method can be terminated at step 340 because the 
25 calibration technique provides direct determination of the absolute values of \x a , and ti' s . 

Also, in other embodiments, the summations in the Equations above need not 
be over every source and detector in the experimental apparatus. For example, source- 
detector measurements between any two sources and any two detectors are sufficient to 
determine all of the source-detector coupling coefficient products corresponding to that 
30 set of sources and detectors. This example would be equivalent to setting N s =2 and N d =2 
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in the above Equations. Once some of the coupling coefficient products are determined 
based on a subset of all possible source-detector measurements, the remaining coupling 
coefficient products can be determined with relatively fewer additional measurements. 

5 Forward and Inverse Calculations 

Techniques for the forward and inverse calculations of light propagation 
within the sample are known in the art. See, e.g., S.R. Arridge in "Optical tomography in 
medical imaging" ibid.. An exemplary formulation of the calculation is described below. 
The absorption coefficient \x a and diffusion coefficients Z>, where 
10 D = u/[3(ja a + jjij )], are expanded into a spatially independent background term \i a ° and 

D 0 , repsectively, and a perturbative spatially dependent terms 5j^a(r) and 5D(r), 
respectively. These terms are then incorporated into the diffusion equation, whose formal 
solution can be expressed as an integral equation by use of the appropriate Green function 
corresponding to the sample's boundary conditions. 
1 5 In the present case, the light energy density is expanded in a perturbative 

series, i.e., U(r) = Uo(r) + U\(r) + and solved to first order. The first-order 
perturbative solution to the heterogeneous equation, in the limit in which U\ is 
given by: 

20 ^ 0 (r J ,r^) = Mexp(/* 0 |r 5 -r rf |)/(47iDo|r J -r,|) (11) 



^1 (r„ r d ) = jj - 5^ (r)uZ)o 1 U 0 (r s , r)G(r, r d ) + V U 0 (r s , r ) - VG(r , r d ) 

yl v o 



d 3 r 



(12) 



25 where M is the amplitude of the source located at r s , G(r,r d ) is the Green function 

solution of the homogeneous equation at detector position r d , u is the speed of light in the 
medium, and ko = [(-ojaa°+ia))/D 0 ] ,/2 is the photon density wave number, where co is the 
source modulation angular frequency. If an infinite medium is assumed, the Green 
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function is G(r y r d ) = exp(/* 0 |r - r^|)/(4n|r - r^|) . The integral in Equation (12) is over 
the entire sample volume. In the embodiment described above, the source is a CW source 
without any modulation. Thus, co goes to zero and ko is purely imaginary, so that both 
U 0 (r) and t/,(r) are real-valued. The first term in the integral of Equation (12) 
5 corresponds to absorbing inhomogeneities, whereas the second term corresponds to 
scattering inhomogeneities. Where it is known that one or other type of inhomogeneity 
dominates, the non-dominant term can be dropped from Equation (12). 

Based on Equations (11) and (12), the forward calculation corresponds to: 

M 



In other embodiments, it is also possible to expand the perturbation series beyond the first 
term. 

For the image reconstruction, i.e., the inverse calculation, the integral in 
15 Equation (12) is digitized into a sum over voxels (i.e., volume elements), and equated to a 

series of the values for U^'Krf/ 1 ), which is extracted from the measurements g(ij) and 
the calibration. This yields a set of coupled linear equations that relates to the values of 
Sji a and SD in each voxel within the sample, which correspond to spatial variations in ^ 
and (is. Many numerical methods are available for solving this system of equations 
20 including, for example, the r algebraic technique Simultaneous Iterative Reconstruction 

Technique (SIRT) and single value decomposition (SVD). For a reference on SIRT, see, 
e.g., A.C. Kak and M. Slaney in Principles of Computerized Tomographic Imaging , 
(IEEE, New York, 1988), Chap. 6, p. 21 1. For a reference on SVD, see, e.g., W.H. Press 
et ah, ibid at Chap. 2, p. 52. 

25 

Software Implementation 

The calibration method can be implemented in hardware or software, or a 
combination of both. The method can be implemented in computer programs using 
standard programming techniques following the method and figures described herein. 
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Program code is applied to input data to perform the functions described herein and 

generate output information. The output information is applied to one or more output 

devices such as a display monitor. 

Each program is preferably implemented in a high level procedural or object 
5 oriented programming language to communicate with a computer system. However, the 

programs can be implemented in assembly or machine language, if desired. In any case, 

the language can be a compiled or interpreted language. 

Each such computer program is preferably stored on a storage medium or 

device (e.g., ROM or magnetic diskette) readable by a general or special purpose 
1 0 programmable computer, for configuring and operating the computer when the storage 

media or device is read by the computer to perform the procedures described herein. The 

computer program can also reside in cache or main memory during program execution. 

The calibration method can also be implemented as a computer-readable storage medium, 

configured with a computer program, where the storage medium so configured causes a 
15 computer to operate in a specific and predefined manner to perform the functions 

described herein. 

For example, referring again to FIG. I, analyzer 140 includes a processor 170, 
and input/output control card 160, a user interface 190 such as a keyboard and monitor, 
and a memory 180. The memory stores a program 185 specifying the steps of the 
20 calibration method. When executed, the program causes the processor to carry out the 
steps of the calibration method. 
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Examples 

Additional aspects, features, and advantages of the invention will now be 
exemplified in the following non-limiting examples. 

5 1 . Diffuse Optical Tomography Apparatus 

A continuous-wave diffuse optical tomography system was built. A schematic 
of the electronic hardware for the system is shown in FIG. 4. 

The optical sources were standard low power diodes, although more powerful 
light-emitting diodes can also be used. The detector array was a monolithic 
1 0 photodiode/preamplifier IC housed in a clear 8-pin DIP package (OPT209 from Burr- 
Brown, Tucson, AZ). Alternative embodiments for the detectors can range from, e.g., 
multianode PMTs to discrete commercial-grade photodiodes with external preamplifiers. 

Optical cross talk was minimized by concealing each detector package in 
opaque heat shrink tubing and providing sufficient separation to further attenuate any 
1 5 stray reflections within the metal housing. Each detector fiber was sheathed in opaque 
tubing as well, which also served to protect the fragile cladding of the polymethyl 
methacrylate (PMMA) fibers from abrasion. Electrical cross talk occurring at the front 
end was not significant since the only high impedance node in the preamp was physically 
removed from the DIP package and the modulation frequency was in the kilohertz range. 
>0 Ground loops were minimized by using an electrically isolated power supply, a battery- 
powered portable laptop computer, and a single-point earth connection to safety ground. 
Cross talk through power lines was minimized with on-board regulation and the use of 
separate power supply decoupling filters at each op-amp. Feedthru among multiple 
channels within the analog multiplexer was minimized by switching only demodulated 
15 (baseband) signals and by reducing the DC currents through the switches by placing a 
buffer directly at the multiplexer output. 

To meet a dynamic range objective of 80dB objective, the system included a 
133 kilosample-per-second, 16-bit A/D converter (Analog Devices AD7884). The high 
sample rate allowed rapid scanning through the 1 6 demodulated outputs during a high 
10 data rate measurements, while still permitting oversampling and subsequent averaging in 
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the digital-domain during slower measurements. The performance goals for the system 
are shown in Table 1 . 



Table 1 

5 ~ 



PARAMETER 


GOAL 


DYNAMIC RANGE 


10,000:1 (80dB) 


NONLINEARITY 


<1% over the 80dB dynamic range 


SETTLING TIME 


<300ms to 0.1% 




^KJ.XJ 1 /o 


DIGITAL RESOLUTION 


1 6 bits 


SOURCE CHANNELS 


9 at 780nm and 9 at 830nm 


SOURCE OPTICAL POWER 


-5mW 


DETECTORS 


16 Si photodiode/preamplifiers (OPT209) 


MODULATION TECHNIQUE 


Single-phase squarewave AM with 
coherent detection 


POSTDETECTION BANDWIDTH 


10to20Hz 


STRAY LIGHT REJECTION 


<1% error under normal illumination levels 


PACKAGING ISSUES 


Must be portable, compact, and extremely 
rugged 


POWER REQUIREMENTS 


120 VAC +/-10%, 50-60Hz 


PATIENT SAFETY 


Leakage current <1jjA, case-to-gnd 
impedance <1Q 



Synchronous detection was used to detect picowatts of source signal under an 
ambient optical background in the microwatt range (as seen by the detector). In 
particular, the laser diodes sources were intensity-modulated with a 50% duty cycle by 

10 modulating the bias current at 2 kHz. Each photodetector preamp output was high-pass 
filtered to remove ambient low frequency signals and then fed into a double-balanced 
mixer. The mixer, which is gated by the same modulator that controls the laser diode 
source intensity, synchronously rectifies the weak modulated source signal to produce a 
DC voltage at the mixer output. All other spurious modulated optical signals (including 

15 those produced by line-powered lamps, computer terminals, multiplexed LED displays, 
etc.) exit the mixer as frequency-shifted AC signals. A subsequent low-pass filter 
strongly attenuates the AC signals, leaving only the DC voltage, which is proportional to 
the magnitude of the detected source energy. The post-detection time constant was set to 
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40 ms and the dwell time was set to about 200 ms to provide sufficient settling for each 
detector channel. A computer was used to select the source and detection channels. 
Table 2 below lists performance results for the system. 

The system noise floor was obtained by calculating the standard deviation of a 
5 number of readings taken at the lowest resolvable source intensity. The noise equivalent 
power was then calculated for an SNR of unity, using the power-per-unit-flux value 
obtained from the linearity measurement. 

Dynamic range and linearity are related: the dynamic range can only be 
defined with reference to a specified linearity limit. The incident beam from an optical 
1 0 source was sampled with a "monitor" fiber that led to a calibrated power meter. The 
linearity was then measured by comparing this value to the system reading over a wide 
range of flux levels. The flux was varied using neutral density filters. A diffuser was used 
to reduce the spatial coherence of the source, which prevented modal noise in the fibers 
from interfering with the measurement. 

15 



Table 2 



PARAMETER 


MEASURED VALUE 


NOISE EQUIVALENT 


<40p W RMS (measured with 32 samples 
per dwell 


DYNAMIC RANGE 


~45,000;1 (92dB) @ 0.75% nonlinearity 
-25,000; 1 (88dB) &). 0.05% nonlinearitv 


LONG-TERM STABILITY 


+/-1% if reading in 30 minutes (half-scale 
output) 


INTERCHANNEL CROSSTALK 


<I:20,000(-86dB) 


STRAY LIGHT REJECTION 


~ DN signal change from darkness to 
normal ambient using cool-white 
fluorescent and incandescent lamps 


TEMPORAL RESPONSE 


~20Hz 


POWER DELIVERED TO TISSUE 


2mW (® 780nm, 8mW (ft 830nm 



Long-term stability was measured by noting the total drift between initial and 
final readings for a fixed probe geometry on a static sample. Interchannel cross talk was 
measured as a single detector channel was alternatively driven between the noise floor 
and full-scale output using the diffuser/ND filter technique described earlier. The largest 
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level change among the other fifteen channels was recorded. Stray light rejection was 
measured by operating the system with a static sample in a dark room and then turning on 
both the fluorescent lights and a computer display located about a meter away. The 
largest change among the sixteen detector channels was recorded. The detected power in 
5 the visible band was ~1 .5 uW. 

The power delivery was measured at the end of -1 meter of 1 mm diameter 
polymethly methacrylate fiber using a calibrated power meter. No index-matching with 
the detector was attempted. This may underestimate the actual power by -4% for in-vivo 
measurements due to index-matching from perspiration trapped at the fiber/skin interface. 

10 

2. Experimental Calibration Results with on Bath of Intralipid and India Ink 

Experimental measurements were performed on a bath of highly scattering 
sample of varying concentrations of Intralipid and India ink using the system described in 
Example 1 . The bath was deep enough to be considered "semi-infinite." The nine 

1 5 sources and sixteen detectors used in the measurement are shown in FIG. 5. The distance 
between neighboring detectors and neighboring sources is 1 cm. The experimental set-up 
was fixed and therefore the coupling coefficients of all channels were nearly constant in 
time although they were different from each other and unknown. 

In a first experiment, the absorption coefficient was incrementally increased 

20 by adding more India ink. The Intralipid concentration, and hence the reduced scattering 
coefficient was kept constant. The calibration method, e.g., steps 300-330 of FIG. 3, was 
then used to determine the-coupling constants, the absorption coefficient for the varying 
concentrations of India Ink, and the constant reduced scattering coefficient. FIG. 6a 
shows the plot of the absorption coefficient determined with the calibration method (filled 

25 squares) as a function of the ink concentration. The results were compared with 

absorption coefficients (open circles) determined without the calibration method, in 
which case the measured signals g(\ j) were fit directly to a homogeneous model for the 
bath. 

Expected results based on independent spectrophotometer measurements are 
30 plotted as a straight line and show the improvement based on the calibration method. For 
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example, at a concentration of 3% of ink solution, the spectrophotometer measurement 
was of 0.052 

cm* 1 , compared to 0.06 cm" 1 and 0.02 cm" 1 , with and without the calibration, respectively. 
The reduced scattering coefficients were determined to be 9.3 cm' 1 and 24 cm' 1 with and 
5 without the calibration, respectively. Once again the result with the calibration was much 
closer to an independent measurement, which was 10 cm" 1 for the reduced scattering 
coefficient. 

Improved accuracy because of the calibration was similarly observed in a 
second set of experiments in which the scattering coefficient was incrementally increased 
10 by adding more India ink. The results for the second set of experiments are shown in 
FIG. 6b. 

3. Simulated Calibration Results 

To determine how the variance in the calculated optical properties depends on 
15 the variance in the coupling coefficients, both with and without use of the calibration 
procedure, simulations were performed to mimic the conditions in Example 2. To 
generate the simulated measurements for the calibration, forward calculations were 
performed with the coupling coefficients of each source and detector randomly varied 
between 0.5 and 1. About 1% Gaussian noise was then added to the forward calculation. 
20 The simulated data was then analyzed with and without the calibration method 

to determine the absorption coefficient and reduced scattering coefficient. The results for 
absorption and constant scattering are shown in the FIGS. 7a and 7b show the results for 
the absorption coefficient and the reduced scattering coefficient, respectively, as 
absorption is being increased, and FIGS. 7c and 7d show the corresponding results as 
25 scattering is being increased. The filled squares indicate results with the calibration, the 
open circles indicate results without the calibration, and the line indicates a fit to expected 
results. The error bars represent the standard deviation caused by varying the coupling 
coefficients in the simulations. The simulations show that the calibration method 
improves both the mean and standard deviation of the determined coefficients. 

30 

4. Simulated Results for a Hidden Object 
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Simulations were performed to apply the calibration method to enhance the 
accuracy of imaging by diffuse optical tomography of an object hidden in a highly 
scattering medium. The background medium was simulated to have a reduced scattering 
coefficient of 10 cm' 1 and an absorption coefficient of 0.05 cm" 1 , and the object was a 2 
5 mm by 2 mm by 2 mm cube centered at a depth of 1 cm under the optical probe array and 
having the same reduced scattering, but double the absorption, as those of the background 
medium. In lateral coordinates, the object was centered along the y-axis and +0.2 cm off- 
center along the x-axis. The positions of the 9 sources and 16 detectors were the same as 
that described in Example 2 (see FIG. 5). The coupling coefficient for each source and 

1 0 detector was randomly varied between 0.5 and 1 and about 1% Gaussian noise was added 
to the simulated measurements. The calibration method was then applied to the simulated 
measurements. The perturbation expansion of Equations (1 1) and (12) and the SIRT 
technique were used to perform the inverse calculation. 

FIGS. 8b and 8a show a reconstructed image of the object with and without, 

1 5 respectively, application of the calibration method to correct the simulated measurement 
prior to the inverse calculation. FIGS. 8a and 8b correspond to one set of coupling 
coefficients and noise. The simulations were repeated with 100 different sets of randomly 
chosen coupling coefficients and noise to observe the variance in the reconstructed 
images. Application of the calibration method reduced the variance by more than 95% in 

20 these simulations. 

5. Experimental Imaging of Hematomas in a Piglet Model 

A piglet model was used in this study. The experimental set-up was a hybrid 
system of diffuse optical tomography and x-ray CT. Measurements were made on the 

25 bench of the x-ray scanner. 

A one-week old piglet, weighted 3 kg, was sedated, intubated, and ventilated 
during the experiment. The femoral artery was catheterized for continuous blood pressure 
monitoring, fluid infusions, and blood extractions. The blood taken from the femoral 
artery was delivered through two small needles inserted 2 cm through scalp, skull, and 

30 brain tissue with a separation of about 2 cm. The combined thickness of the scalp and 

skull is about 1 cm, so the blood was injected into the brain at a depth of about 1 cm. The 
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injection speed was controlled by a step motor at a rate of 42 u,l/minute. Total 625uJ 
blood was injected for each bleed in 15 minutes respectively. Images were reconstructed 
every 15 seconds. Bleed A was created first, then bleed B. The optical probe was placed 
on top of the piglet's head. 
5 Referring to FIG. 9, the probe has 16 detectors (large circles) and 9 sources at 

each of 780 nm and 830 nm. The positions of bleeds A and B relative to the probe are 
also shown in FIG. 9. Before starting to inject the blood into A, a baseline was measured 
and used to find the coupling coefficient of each source-detector channel and also the 
background optical properties. The calculated absorption and effective scattering 

1 0 coefficients were 0.0672 cm" 1 and 8.44 cm" 1 , respectively, at 780 nm, and 0.0666 cm"' and 
7.61 cm 1 , respectively, at 830 nm. These values correspond to a S0 2 of 58% and HbT of 
78 jil/mol. After bleed A was generated, another baseline was measured before injecting 
blood into B. X-ray CT images were taken after the optical measurements to identify the 
positions of the bleeds. 

1 5 2D images were reconstructed based on the DOT measurements using 

simultaneous iterative reconstruction technique (SIRT) and assuming that the piglet head 
is semi-infinite. FIGS. lOa-d shows the time-course of the reconstructed images of bleed 
A at both wavelengths, both with and without using the calibration method to correct the 
DOT measurements. FIGS. lOa-d show that the application of the calibration method 

20 greatly reduces the presence of artifacts in the bottom and right sides of the images. The 
calibration method also improves the image amplitude. 

FIGS. 1 la-b show a time-course plot of the peak image value of bleed A at the 
two wavelengths. Without calibration (FIG. 1 la), the plot is noisy with no clear build up 
of the bleed. After the calibration (FIG. 1 lb), however, the development of the bleed is 

25 clearly visible. 

FIG. 12 shows time-course images of both bleeds A and B as the volume of 
bleed B was increased. The base line and coupling coefficient were the same as those 
used in reconstructing images of bleed A. The respective intensities of the A and B 
bleeds differ because their depths differ. 
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Other Embodiments 
It is to be understood that while the invention has been described in 
conjunction with the detailed description thereof, that the foregoing description is 
intended to illustrate and not limit the scope of the invention, which is defined by the 
5 scope of the appended claims. 

Other aspects, advantages, and modifications are within the scope of the 
following claims. 

What is claimed is: 

10 
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1 . A system for making optical measurements of a sample, the system 
comprising: 

at least two optical sources which during operation couple optical radiation 
into the sample at spatially separated locations; 
5 at- least two optical detectors positioned to receive optical radiation emitted 

from the sample at spatially separated locations in response to the optical radiation from 
the sources, wherein the signal g(ij) produced by the /» det e C tor in response to the 
optical radiation from the source can be expressed as g(i,j) = S'D'/fcy), where 
f(ij) depends only on the properties of the sample, ^ is a coupling coefficient for the fi 
1 0 source, and D) is a coupling coefficient for the /h detector; and 

an analyzer which during operation calculates the value of the product D k , 
wherein 5> is a coupling coefficient for the 1 th source and D k is a coupling coefficient for 
the k* detector, for at least one of the source-detector pairs denoted by the superscripts I 
and k based on the signals produced by the detectors and simulated values of f(ij) 
1 5 corresponding to a model of the optical properties of the sample. 

2. The system of claim 1, wherein during operation the analyzer calculates the 
value of the product S«Z» for every source-detector pair based on the detector signals and 
the simulated values of f(ij). 



20 



3. The system of claim 2, wherein the analyzer further calculates experimental 
values oifdj) based on the calculated values of &DS and the signals g(ij) using the 
expression g(i,j) = S'D^f(i,j), and performs an inverse calculation on the experimental 
values for f(ij) to determine spatial variations in at least one optical property of the 



25 sample. 



4. The system of claim 3, wherein the at least one optical property comprises 
at least one of an absorption coefficient and a reduced scattering coefficient. 
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5. The system of claim 3, wherein during operation the analyzer modifies the 
model of the sample based on the determined spatial variations and repeats the calculation 
of the values of the product SPD* for every source-detector pair using the modified model. 

5 6." The system of claim 1, wherein the model corresponds to the sample being 

homogeneous. 

7. The system of claim 1, wherein the analyzer simulates the values oif(iJ) 
according to the expression: 



10 



fO'J) ^exp^-to Jl^lj 



where Jr^-j is the distance between the i th source and the j lh detector, where jj' s is the 
reduced scattering coefficient, and where /u a is the absorption coefficient. 

1 5 8. The system of claim 1, wherein the analyzer calculates the value of the 

product S'D^ by minimizing the expression: 



where 



WJ,k,l) = A^Af 



20 A , k _ N s g g (ij) 

A ji = ^y gOJl 



and N s is the number of sources and N D is the number of detectors. 
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9. The system of claim 8, wherein the model corresponds to the sample being 
homogeneous, and wherein during operation the analyzer calculates at least one of the 
absorption coefficient \x a and the reduced scattering coefficient \x' s by,minimizing 

5 F(s l D k ) with respect to the product S l D k and the at least one of the absorption and 

scattering coefficients, f(s 1 D*) implicitly depending on \x a and \x! s through /ft/J. 

10. The system of claim 9, wherein during operation the analyzer calculates 
both of the absorption and scattering coefficients by minimizing f{s 1 D k ) with respect to 

1 0 the product sLD k and the absorption and scattering coefficients. 

1 1 . The system of claim 8, wherein the analyzer calculates the product S^D™ 
for every source-detector pair by minimizing F{s m D n \ 

15 1 2 - The system of claim 8, where the analyzer calculates the product S*D) for 

every remaining source-detector pair according to S'D j = ^LlihiH 

S'D k ' 

13. The system of claim 1, wherein g(ij),f(ij), 5% and D) are all real-valued. 

20 14- The system of claim 1, wherein the sources couple continuous- wave 

optical radiation into the sample. 

15. A method for calibrating an optical measurement system including at least 
two optical sources and at least two optical detectors, wherein the sources couple optical 
25 radiation into a sample at spatially separated locations and the detectors are positioned to 
receive optical radiation emitted from the sample at spatially separated locations and 
generate signals in response to the optical radiation from the sources, the method 
comprising: 



-30- 



WO 01/19241 



PCT/US00/25467 



providing the signals generated by the detectors, wherein the signal g(ij) 
generated by the fi 1 detector in response to the optical radiation from the /* source can be 
expressed as g(i,j) = S* D J f(i 9 j), where f(ij) depends only on the properties of the 
sample, 5* is_a coupling coefficient for the source, and £U is a coupling coefficient for 
5 the yth detector; and 

calculating the value of the product S'z^ for at least one of the source-detector 
pairs based on the signals generated by the detectors and simulated values of f(ij) 
corresponding to a model of the optical properties of the sample. 

10 16. A computer readable medium comprising a program which causes a 

processor to perform the steps of claim 15. 

1 7. The method of claim 1 5, wherein the value of the product D k is 

calculated for the product S*DS for every source-detector pair based on the detector 
15 signals and the simulated values of f(ij). 

18. The method of claim 17, further comprising calculating experimental 

values of f(ij) based on the calculated values of &DS and the signals g(ij) using the 

expression g(i,j) = S'D-? f(i,j), and performing an inverse calculation on the 

20 experimental values for f(ij) to determine spatial variations in at least one optical 
property of the sample. 

19. The method of claim 15, wherein values of f(ij) are simulated according 
to the expression: 
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where \r, y \ is the distance between the i* source and the j* detector, where M ' M is the 
reduced scattering coefficient, and where Mo is the absorption coefficient. 
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